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Hybrid density functionals show great promise for chemically-accurate first principles calculations, 
but their high computational cost limits their application in non-trivial studies, such as exploration 
of reaction pathways of adsorbents on periodic surfaces. One factor responsible for their increased 
cost is the dense Brillouin-zone sampling necessary to accurately resolve an integrable singularity 
in the exact exchange energy. We analyze this singularity within an intuitive formalism based on 
Wannier-function localization and analytically prove Wigner-Seitz truncation to be the ideal method 
for regularizing the Coulomb potential in the exchange kernel. We show that this method is limited 
only by Brillouin-zone discretization errors in the Kohn-Sham orbitals, and hence converges the 
exchange energy exponentially with the number of fc-points used to sample the Brillouin zone for all 
but zero-temperature metallic systems. To facilitate the implementation of this method, we develop 
a general construction for the plane-wave Coulomb kernel truncated on the Wigner-Seitz cell in 
one, two or three lattice directions. We compare several regularization methods for the exchange 
kernel in a variety of real systems including low-symmetry crystals and low-dimensional materials. 
We find that our Wigner-Seitz truncation systematically yields the best fc-point convergence for the 
exchange energy of all these systems and delivers an accuracy to hybrid functionals comparable to 
semi-local and screened-exchange functionals at identical fc-point sets. 



I. INTRODUCTION 

Density-functional theory-^ forms the basis for ab 
initio theoretical studies of the ground-state electronic 
structure of materials, and serves as the starting point 
for many-body perturbation theories such as GW^ and 
Bethe-Salpeter equatioiP (BSE) calculations. Standard 
semi-local approximations to exchange and correlation in 
density-functional theory, such as the local-density and 
generalized-gradient approximations, are remarkably ac- 
curate for a variety of properties such as lattice constants, 
equilibrium geometries and elastic moduli, but are not 
sufficiently accurate for the energetics and kinetics of 
chemical reactions.^ 

Hybrid density functionals, which replace a fraction of 
the approximate semi-local exchange energy with the ex- 
act non-local Fock exchange energy improve upon the 
accuracy of semi-local functionals and have been widely 
applied for first-principles thermo-chemistryP Variants of 
these functional^ enable calculations for solids and sur- 
faces with accuracy sufficient for predicting atomic-scale 
processes at room temperature. However, hybrid func- 
tionals require a greater number of fc-points than semi- 
local functionals for comparable accuracy in Brillouin- 
zone discretization for periodic systems. This increases 
their already high computational cost and limits their ap- 
plicability, so that practical studies of surface reactions 
and phenomena such as catalysis remain tantalizingly out 
of reach. 

The need for finer fc-point sampling in hybrid func- 
tionals stems from the Brillouin-zone integrals over the 
singular Coulomb kernel in the exact exchange energy of 
periodic systems. Similar operators also appear in GW 
and BSE calculations, and so successful methods to ad- 
dress this issue have implications for excited state meth- 



ods as well. 

The discretization error in singular reciprocal space in- 
tegrals critically depends on the technique used to han- 
dle the singular contributions. The standard auxiliary- 
function approacfP^ replaces the divergent terms with 
an average around the singularity computed using an 
auxiliary function with the same singularity as the 
Coulomb kernel. However, such methods that replace 
only the divergent terms lead to polynomial convergence 
(N^ 1 ) with the number of fc-points iV^ES compared to 
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the standard exponential convergence (exp(— aN k )) in 
the case of semi-local functionals. 

Some hybrid functionals 11 achieve fc-point convergence 
comparable to fully semi-local functionals by replacing 
the exact exchange energy with a screened exchange en- 
ergy computed using a short-ranged kernel erfc(o;r)/r 
that is not singular at long wave- vectors. The remaining 
long-ranged part erf (ujr)/r is treated using a semi-local 
(generalized-gradient) approximation. The predictions of 
these functionals depend on the screening parameter co 
due to this additional approximation, and are less accu- 
rate than those of exact-exchange hybrids for some prop- 
erties such as elastic constants of periodic systems!^ A 
method with comparable computational cost but with 
the exact 1/r kernel would therefore be highly valuable. 

One approach that shows promising results for the ex- 
act exchange energy avoids the singularity by imposing 
a real-space cutoff on the Coulomb kernel with a length- 
scale dependent on the fc-point meshP^HIH The spherical 
truncation employed in this approach, however, works 
well only for high-symmetry crystals, whose Wigner-Seitz 
cells more or less resemble spheres. The valuable gains 
afforded by such an approach in these special cases indi- 
cates the need for a detailed understanding of why trun- 
cated potentials work for calculating the exchange energy 
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of periodic systems so as to point the way to a general 
method applicable to all systems. 



Following this program, Section II A analyzes the sin- 
gularity in the exact exchange energy of periodic sys- 
tems using a formalism based on Wannier-function lo- 
calization. This allows us to prove analytically that 
Wigner-Seitz truncation of the Coulomb potential is 
the ideal regularization method with accuracy limited 
only by Brillouin-zone discretization errors in the Kohn- 
Sham orbitals themselves; the appendix develops a set 
of techniques necessary to truncate the Coulomb po- 
tential on Wigner-Seitz cells. Section |II B| then gen- 
eralizes regularization methods for the exchange en- 
ergy in three-dimensional systems to slab-like (two- 
dimensional) and wire-like (one-dimensional) geometries 
using partially-truncated potentials. Finally, Section |III| 
compares truncated-potential and auxiliary-function ap- 
proaches for a variety of real materials with high- and 
low-symmetry crystal systems, electronic structure rang- 
ing from insulating to metallic, and dimensionality rang- 
ing from three to one. 

The results indicate that employing Wigner-Seitz trun- 
cated Coulomb kernels is systematically more accurate 
than other methods for dealing with the singularity in 
Fock exchange. Moreover, at equivalent Brillouin-zone 
sampling, the accuracy of hybrid functionals which in- 
clude exact exchange computed using the Wigner-Seitz 
truncated method rivals that of functionals that only 
include screened exchange, and even that of semi-local 
functionals which include no non-local exchange contri- 
butions whatsoever. 



II. EXCHANGE IN PERIODIC SYSTEMS 

The exact Fock exchange energy of a finite system with 
Kohn-Sham orbitals ipi a (r) and corresponding occupa- 
tion numbers fi a is given by the non-singular expression, 



(1) 



Here, we work with atomic units e 2 /(47reo) = ti 2 /m e = 1, 
so that the unit of distance is the Bohr (ao) and the unit 
of energy is the Hartree (Eh). The exchange energy is 
always a sum of independent contributions for each spin 
channel; the rest of this section omits the spin index a 
and the implicit sum over a for clarity. 

In a periodic system, the Kohn-Sham orbitals take the 
Bloch form ipf(r) = e r u k (r), where u k (r) are periodic 
functions normalized on the unit cell of volume f2 and 
labeled by band index i as well as a wave-vector fc in the 
Brillouin zone. The exchange energy of the periodic 



system per unit cell (per spin) is 
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where J BZ denotes integration over the Brillouin zone. 
The conventional treatment of this energy begins with a 
plane-wave expansion of the product densities 
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where 5 kk ~ are the Fourier components of those densi- 
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ties at reciprocal lattice vectors G. This treatment then 
rewrites ^ in Fourier space as 
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where the periodic Coulomb kernel = Air/q 2 , so that 
for each k' at G = 0, the integral over q = k — k' is 
singular at q = 0. This singularity is integrable since 
near q = 0, the integral ~ J Airq 2 dq^ I . 

The above approach, however, is problematic for any 
practical calculation where Brillouin-zone integrals are 
approximated using a finite quadrature, that is, as a 
weighted sum over a set of 'fc-points'. In this paper, we 
restrict our attention to the commonly employed Gauss- 
Fourier quadratures, which correspond to uniform fc- 
point meshes such as the Monkhorst-Pack gridPSl The 
exchange energy computed in practice is therefore 
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where is the total number of fc-points used for Bril- 
louin zone sampling. In principle, the exchange energy 
would converge with increasing density of fc-points even 
if the singular terms are dropped, or equivalently, the 
Coulomb kernel is regularized with K q= o — 0, as usual. 
However, that results in an 0(5q) error, where Sq is 
the typical distance between neighboring fc-points, which 
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leads to an impractically slow N k convergence. 

Auxiliary-function methods 8 address this poor Bril- 
louin zone convergence of the Fock exchange energy by 
choosing a value for the G — 0, k = k' term in (j5j) 
that captures the average contribution of 47r/|fc' — fc| 2 
in the neighborhood of k — fc'. These methods correct 
for the finite quadrature error by setting this term to 
the difference between the exact integral and the dis- 
crete fc-point sum over the Brillouin zone of a function 
/(g) that matches the periodicity and the An/q 2 singular- 
ity of the integrand in Q. For a uniform fc-point mesh, 
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this amounts to replacing the Coulomb kernel Kq in (|5 
with 



K 9 



= Uir/q 2 , q^O 



(6) 



where the discrete sum runs over 5k in the fc-point dif- 
ference mesh excluding the T point 

The original method of Gygi and Baldereschi presented 
such a function for the face-centered cubic lattice which 
could be integrated analytically. Wenzien and coworkers 9 
constructed similar functions for a few other lattice sys- 
tems and tabulated the corresponding q = corrections 
numerically. Carrier and coworkers^ constructed a gen- 
eral function that works for all Bravais lattices and pre- 
scribed a general scheme for computing the Brillouin zone 
integral. Below, when making comparisons to our trun- 
cated potential method, we employ this last variant of the 
the auxiliary-function method due to its generality, and 
refer the reader interested in further details of auxiliary- 
function methods to Ref. [TO] 

With the correctly chosen G = term for k = k' , the 
auxiliary function methods achieve N^ 1 convergence^ in 
the exchange energy. Duchemin and GygP^have general- 
ized the method to achieve Nj7 2 convergence by introduc- 
ing 'curvature corrections' which, in our notation above, 
amount to correcting the Coulomb kernel at q = as well 
as q that correspond to nearest neighbor displacements 
in the /c-point mesh. For many systems, this method 
yields reasonable accuracy with modest /c-point meshes; 
however the asymptotic polynomial convergence of the 
exchange energy is still slower than the exponential con- 
vergence of the total energy of semi-local density func- 
tionals. We next describe a method for the exchange 
energy that achieves this exponential convergence. 



A. Real-space analysis of asymptotic convergence 

An alternate scheme to improve the Brillouin zone con- 
vergence of exact exchange imposes a real-space cutoff 
on the Coulomb kernel in ([5]). This scheme been shown 
to work reasonably well for high-symmetry crystals,^ 
but the reasons for its success remain somewhat myste- 
rious. Two possible explanations have been offered. The 
first is that the method satisfies the normalization con- 
straints of the exchange hole by appropriately truncating 
the Coulomb potential.^ The second is that the effective 
distinguishability of electrons amongst different fc-point 
sampled supercells requires suppression of the exchange 
interaction between supercells!^ These explanations do 
not elucidate the underlying reason for an infinite-range 
interaction to be best numerically approximated by a 
finite-range one nor specify what form that finite-ranged 
interaction should take, and they do not lend themselves 
to an analysis of the accuracy or convergence properties 
of such an approximation. To provide such an explana- 
tion and to identify the ideal form of the truncation, we 



now analyze the Fock exchange interaction computed on 
finite fc-point meshes in real space, and show that the 
need for truncating the Coulomb potential arises both 
naturally and in a particular form. 

We start by rearranging the exchange energy of the 
periodic system ^ as 
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a sum of Coulomb self-energies of the pair densities 
pf-{r) = (r)w^(r) of the Wannier-like functions 
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Indeed, order- N calculations of the exchange energy^ us- 
ing maximally-localized Wannier functions 20 employ sim- 
ilar transformations. In contrast, we use Q, which is ex- 
actly equivalent to ([2| , only as a tool to analyze standard 
reciprocal-space methods. 

In the case of insulators, where /* = 1 for all occupied 
bands, the Wannier-like functions of ^ are just the stan- 
dard Wannier functions and thus are exponentially 
localized around the sites RW^ For the case of metals, 
where the occupations /* are not constant, the wf are 
linear combinations of Wannier functions localized on dif- 
ferent lattice sites 
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with coefficients Fi(R' — R) given by a Fourier transform 
of the square-root of the band occupation. For metals at 
zero temperature, the occupations are discontinuous at 
a Fermi surface which leads to a polynomially decaying 
F t . In particular, for metals in three dimensions with 
a compact two-dimensional Fermi surface, Fi(R) ~ Rr 2 
for large R. At finite temperature T, Fi also decays ex- 
ponentially with a length scale inversely proportional to 
T. Consequently, the wf in metals are localized polyno- 
mially at T = and exponentially with a temperature- 
dependent decay length at T ^ 0. In all these cases, the 
localization of closely mirrors that of the one-particle 
density matrix.^ 

Given the above properties of the Wannier-like func- 
tions wf-, each pair density pfj{r) is localized and di- 
minishes with increasing R, since it is a product of one 
function localized around and another around R. The 
magnitude of the pair density decreases exponentially for 
insulators and as R~ 2 for metals at zero temperature, so 
that the corresponding Coulomb self-energies decay ex- 
ponentially and as R~ respectively. The sum over unit 
cells in (JtJ) converges in all these cases. 
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Next, to understand the convergence properties of ac- 
tual calculations, we repeat the above transformations 
with a finite fc-point mesh with samples instead of 
continuous integrals over the Brillouin zone. A uni- 
form fc-point mesh centered on the T-point corresponds 
to Kohn-Sham orbitals that are periodic on a supercell 
of volume N^n. For uniform off-r meshes, such as the 
Monkhorst-Pack grid, the orbitals share a common Bloch 
phase on an iVj.fi supercell. In all these cases, the pair 
densities pfj{r) are periodic on that Nf.0, supercell. Be- 
cause of this periodicity, the expression for the exchange 
energy Q then remains unmodified except that one of 
the integrals over space is restricted to a single iVfcfi su- 
percell. Converting the other integral over all of space 
to a sum over integrals restricted to each iVfcfi supercell, 
and using the periodicity of the pair density, 

Ex = ^J2 [ W d?pf 3 *{r)K(r- 
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where S are lattice vectors of the effective superlat- 
tice of cell volume iVfcfi which arises from finite fc- 
point sampling. The summand in (11) falls off only 
as 1/S causing K(r) to diverge for all 

K(r) — 1 V- -* p i(k+G)-f_ 



In fact, 
so that it is again the 



'k,G" \k+G\ 2 

k = G = component which needs special handling, as 
above. 

Now consider a sufficiently dense fc-point mesh such 
that the corresponding supercell is much larger than the 
spatial extent of the localized wf. In that case, the peri- 
odic versions of pfj at finite fc-point sampling are identical 
(with exponentially small errors) in one, appropriately 
centered, supercell to the original non-periodic localized 
ones. Therefore, the contribution from the 5 = term 



in K(r) to ( 10 ), apart from errors decaying exponentially 
with the density of the fc-point mesh, is the true exchange 
energy of the infinite system Q. The non-exponentially 
decaying errors in ( 10 1 arise from the contributions due 



to all the other super-cells S ^ 0, which thus can be elim- 
inated completely by truncating the Coulomb potential 
so that K(r) = l/r\ 

Such truncation of the Coulomb potential on the 
Wigner-Seitz cell of the fc-point sampled superlattice with 
supercell volume iVfcfi, in practice simply amounts to re- 
placing Kq in the standard reciprocal-space expression 
([5]) by the Fourier transform of the truncated potential. 
The minimum image convention (MIC) algorithm^ em- 
ploys 
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FIG. 1. Comparison of the discrepancy in the effective 
Coulomb Kernel from 1/r of various methods for comput- 
ing the Fock exchange energy. This discrepancy is plotted 
along the x-direction from the origin to the boundary of the 
cubic Wigner-Seitz cell of the fc-point sampled super-lattice 
of a simple cubic lattice of lattice constant a. Both axes have 
been scaled to be dimensionless. The truncated Coulomb po- 
tentials have zero discrepancy for most or all of the super- 
cell. The auxiliary-function methods shift the periodic kernel 
to minimize the overall discrepancy, while the probe-charge 
Ewald-sum method pins the discrepancy at the origin to 0. 



and enables efficient construction of truncated kernels in 
the plane- wave basis. The appendix derives this algo- 
rithm for arbitrary lattice systems as equation ( |A4[ ); see 
appendix |A 1| for a detailed explanation of all the terms 
and approximations involved in (12 1. When such a trun- 



(12) 



cated kernel replaces the periodic Coulomb kernel, the 
remaining error in the exchange energy is due to the de- 
viations of the periodic pfj from the infinite ones, within 
one supercell, as discussed above. These deviations de- 
cay exponentially when the wf' are exponentially local- 
ized, leading to exponential convergence of the exchange 
energy with the number of fc-points. 

The general analysis we developed above not only al- 
lows us to establish the Wigner-Seitz super-cell trun- 
cated potential as the natural method with ideal asymp- 
totic convergence, but also provides the framework for 
establishing analytically the convergence properties of 
other methods a priori by simply comparing the effec- 
tive Coulomb kernels which they employ. Ultimately, 
each reciprocal-space method prescribes a kernel K(q) 
for use in ([5]), which translates to some K(r) in the real 
space version ( |10| . The error in any method compared 
to the ideal case is governed by the discrepancy of K{r) 
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from 1/r, given by 
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(13) 



for a particular fc-point set. 

FIG. [I] compares the discrepancies in the effective 
Coulomb kernels of various methods for a cubic lattice 
over a radial slice of the effective supercell with different 
/c-point meshes. 6K is exactly zero within the domain 
of truncation for any truncated Coulomb potential, and 
hence it is zero for the Wigner-Seitz truncated kernel in 
the entire supercell. Truncation on a sphere of volume 
iVfcfi, as proposed by Spencer and Alavip^ achieves expo- 
nential convergence as well, but exhibits its asymptotic 
properties starting at larger Nk since the sphere does not 
tile with the super-lattice and leads to some overlap be- 
tween supercells. 

The standard periodic kernel with G = projected 
out has a large discrepancy ~ 1/L } where L is the linear 
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dimension of the supercell, and hence exhibits ~ N k 
convergence. The auxiliary-function methods add a con- 
stant offset to the periodic kernel by adjusting the G = 
component, which is chosen to minimize the discrepancy 
over the supercell with some weight which depends on 
the choice of auxiliary function f(q). The curvature cor- 
rections proposed by Duchemin and Gygf^ adjust the 
constant offset in addition to the coefficients of e lSk ' r for 
the nearest-neighbors 6k in the fe-point mesh. For the 
cubic supercell of length L plotted in FIG. [T] this cor- 
rection happens to be (cos(27rx/L) — l)/(nL) in the x 
slice. This additional freedom reduces the magnitude of 
the discrepancy more rapidly with increasing Nk than 
the case when only G = is adjusted. 

An interesting alternative to the auxiliary-function 
method with similar accuracy is the probe-charge Ewald- 
sum method.^ Here, the G = component of the 
Coulomb kernel is set to the potential at the origin from 
an array of unit negative charges placed at all points 
of the fc-point sampled super-lattice except the origin, 
a Coulomb sum which can be computed readily using 
the Ewald methodP^ This amounts to adjusting the dis- 
crepancy at the origin, 6K(0), to zero, as can be seen 
in FIG. [I] Alternately, this method attempts to cancel 
the contributions from 5^0 supercells of (111 in (10) 



by neutralizing all those supercells with a point charge. 
This would be exact if each p was spherically symmetric, 
but in reality incurs an error asymptotically dominated 
by dipole-dipole interactions ~ L~ 3 and hence converges 
as NZ with fc-point sampling. 

The Wannier-function formalism presented here clearly 
establishes the advantage of the truncated-potential 
methods and elucidates the asymptotic convergence of 
the exchange energy computed using different methods. 
Note that none of these methods require Wannier func- 
tions in practice; each method prescribes a different re- 
placement for the periodic plane-wave Coulomb kernel 
K q = Air/q 2 in the standard reciprocal-space expression 



(§. We compare the accuracy of these methods and 
demonstrate their analytically-predicted asymptotic be- 
havior for real materials in Section Hill 



B. Extension to low-dimensional systems 

The preceding section shows how the integrable sin- 
gularity in reciprocal-space calculations of the exchange 
energy of systems with three-dimensional periodicity can 
be regularized using auxiliary function methods or, ide- 
ally, by truncating the Coulomb potential to the Wigner- 
Seitz cell simply by modifying the Fourier transform of 
the Coulomb kernel. Reciprocal-space methods can also 
be applied to systems with lower-dimensional periodic- 
ity by truncating the Coulomb potential along a subset 
of lattice directions. The appendix details these types 
of truncations as well. Specifically, the exchange energy 
in these geometries can be computed using ([5| by em- 
ploying a part ially truncated Coulomb potential given by 



(A6) or (A8) for Kg; and restricting the Appoints sums 



to the two- or one-dimensional Brillouin zone of the peri- 
odic directions alone. The exchange energy still contains 
an integrable singularity, q^ 1 for slabs and lnq for wires, 
that again needs to be addressed just as it did for bulk 
systems. Before going on to present results in the next 
section, here, we briefly describe how we extend each of 
the methods discussed above for bulk systems to the cases 
of these lower-dimensional geometries. 

First, the auxiliary-function method for these geome- 
tries replaces the G = component of the appropriate 
partially truncated Coulomb kernel to obtain 
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where fin is the area/length of the two-/one-dimensional 
unit cell in the periodic directions, BZy is the correspond- 
ing Brillouin zone and is the length/area of the ar- 
tificial periodicity along the truncated directions. The 
auxiliary function needs to match the singularity of the 
truncated Coulomb kernel, and we adapt Carrier and 
coworkers' construction for arbitrary three-dimensional 
lattices^ to lower dimensions. For a slab with lattice 
basis vectors a\ and cii in the periodic directions, and 
corresponding reciprocal lattice vectors b\ and &2, the 
function 
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is periodic on the reciprocal lattice with a single singu- 
larity in the Brillouin zone at q = for arbitrary lattice 
vectors. Similarly, 



f(q) = -2 7 + In - 
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(16) 
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is a suitable auxiliary function for a wire with lattice vec- 
tor a along the periodic direction, where 7 is the Euler- 
Mascheroni constant. In this case, both the integral and 
the sum in ( 14 ) can be performed analytically to yield 
K^ = 2VL ± {\og2N k a- 1 ). 

Next, the probe-charge Ewald compensation method 
generalizes trivially to the slab and wire geometries, 
and only requires the substitution of the usual three- 
dimensional Ewald sum with the appropriate lower- 



dimensional analog. (See Appendix A3 for details.) In- 
terestingly, this method yields the same G — compo- 
nent for the wire-geometry exchange kernel as the auxil- 
iary function method above, so that the two methods are 
identical for wires. 

Finally, the Wigner-Seitz supercell truncated Coulomb 
potential requires no modification for partially periodic 
systems. The truncation domain remains the Wigner- 
Seitz cell of the fc-point sampled super-lattice; one or 
two lattice directions have only a single fc-point and the 
boundaries of the supercell coincide with the unit cell 
in those directions. These Wigner-Seitz cells become in- 
creasingly anisotropic with increasing Nk and spherical 
truncation is no longer a viable option. 



Method 


Compute Time [s] 


Wigner-Seitz truncated 


555 ± 11 


Spherical truncated 


874 ± 12 


Auxiliary function 


552 ± 11 


Probe-charge Ewald 


543 ± 10 



TABLE I. Comparison of the average computation time for 
the exact exchange energy using different regularization meth- 
ods for hexagonal silicon carbide with 8x8x8 fc-point sam- 
pling. The timing statistics are from ten calculations for each 
method on identical 12-core Xeon compute nodes. 



change energy puts the convergence of hybrid function- 
als employing the exact non-local exchange energy (e.g. 
PBEG 7 ) on par with that of the screened-exchange func- 
tionals (e.g. HSE06) and even that of semi-local func- 
tionals employing no non-local exchange whatsoever (e.g. 

pbeP). 



A. Computational Details 



III. RESULTS 

The analysis of Section [TT] identifies Coulomb trunca- 
tion with its asymptotic exponential convergence as the 
natural choice for computing the exchange energy of pe- 
riodic systems, in comparison to auxiliary-function meth- 
ods with polynomial convergence. Here, we compare 
the accuracy of all these methods and demonstrate their 
analytically-established asymptotic behavior for real ma- 
terials with a variety of electronic structures and dimen- 
sionalities. 

Specifically, we consider four methods for computing 
the exact exchange energy, the Wigner-Seitz truncated 
potential introduced here, the spherical truncation of 
Spencer and Alavij^ the probe-charge Ewald compensa- 
tion method,^ and the auxiliary-function method with 
the general function applicable to all lattice systems by 
Carrier and coworkers The first three of these meth- 
ods trivially generalize to lower dimensions, while the 
auxiliary-function method requires minor modifications 
as detailed in Section Hi Bl 

We also compare the convergence of the exact ex- 
change energy using the above methods to that of the 
erf-scrccned exchange employed in the range-separated 
HSE06 hybrid functionalP3 In this functional, the 
Coulomb kernel in the non-local exchange energy is re- 
placed by the short-ranged erfc(wr)/r with lj = O.lla^" 1 , 
while the long-ranged part is approximated using a semi- 
local functional. The screened exchange avoids the G = 
singularity and the HSE06 functional has so far achieved 
superior fc-point convergence compared to regular hybrid 
functionals with exact exchange.^ Here, we demonstrate 
that employing Wigner-Seitz truncation for the exact ex- 



We have implemented all these methods in the open- 
source plane-wave density-functional software JDFTxf22l 
where they are now publicly available. The specifics 
of these implementations are that the auxiliary-function 
and probe-charge Ewald methods simply replace only the 
G = value of the k — kl Coulomb kernel in ^ with 
a precomputed value. The truncated potential methods, 



on the other hand, alter K$ in ([5J) for all q = G + k! — k. 
Spherical truncation uses Kq defined analytically via 
( A2 ) . Wigner-Seitz truncation employs a precomputed 



kernel calculated by applying the MIC algorithm (A4) 
on the supercell, as detailed in the appendix, and then 
redistributing the resulting supercell kernel to unit cell 
kernels for each k! — k. 

TABLE [I] shows that the computational overhead for 
looking up the precomputed kernel in the Wigner-Seitz 
truncated method is negligible, and results in compute 
times equal to the auxiliary function and probe-charge 
Ewald methods, within run-to-run variations. In fact, 
this overhead is negligible compared to that of the ex- 
tra transcendental (cosine) evaluations in spherical trun- 
cation; precomputing the kernel also optimizes spheri- 
cal truncation and we report the analytical evaluation 
time here only to illustrate the negligible lookup over- 
head. Next, the computational effort to calculate and 
predistribute the kernel is negligible compared to a sin- 
gle evaluation of the exchange energy: a mere 1.4 s for 
the example of TABLE [T] Finally, the memory overhead 
of the precomputed kernel is comparable to four Kohn- 
Sham bands, and is therefore negligible for most systems. 

In order to study a large number of materials and fc- 
point configurations within the available computational 
resources, rather than performing fully self-consistent 
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System 


Unit cell 


N a 


a [A] 


c[k] 


Ref. 


FIG. 


2H-SiC 


Hexagonal 


4 


3.076 


5.048 


EH 


3 9 


'b) 


3C-SiC 


FCC 


2 


4.3596 


_ 


EH 


l_ - 


b) 

=1 


Ice XIc 


BCT a 


6 


4.385 


6.219 


E3 
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Si 


FCC 


2 


5.431 


_ 


E3 [2 




|>) 


Platinum 


FCC 


1 


3.924 


- 


M 


6 9 




Diamond 


FCC 


2 


3.567 
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f 


Graphite 


Hexagonal 


4 


2.461 


6.709 


Ml 


[j 




Graphene 


Hexagonal 


2 


2.46 


10 b 


PBE' 






(8,0) SWCNT 1 * 


Tetragonal 32 


25^ 


4.32 


PBE C 


r 





122 



a Body-centered tetragonal. 

6 Coulomb potential is truncated along these directions. 

c DFT lattice constants using the PBE functional.^ 

d Single-walled carbon nanotube. 

TABLE II. Unit-cell parameters for the systems studied here, 
including citations for experimental lattice constants and ref- 
erences to figures with corresponding results. N a is the num- 
ber of atoms in the primitive unit cell of each calculation. 



calculations, we first determine converged Kohn-Sham 
orbitals of a density-functional calculation using the 
semi-local PBE exchange and correlation functional,^ 
and then compute the exchange energies from these or- 
bitals according to the above methods. The calculations 
employ norm-conserving pseudopotentials at a kinetic en- 
ergy cutoff of 30 Eh- TABLE |n] summarizes the unit-cell 
parameters for the systems studied below. 

Figures [2][8] show the deviation of the exact and 
screened exchange energies at finite fc-point configura- 
tions from their fc-point-converged values for a variety of 
systems. The left-hand panels show this deviation for 
coarse fc-point meshes on a linear scale, while the right- 
hand panels illustrate the asymptotic convergence on a 
logarithmic energy scale. The base line of fc-points is in- 
sufficient to reliably fit power laws and the dotted lines 
are only a visual guide with the expected exponent for 
polynomial convergence. For the three-dimensional sys- 
tems, we study both isotropic and anisotropic fc-point 
meshes. The plots explicitly label anisotropic fc-point 
configurations, whereas the unlabeled points at integer 



values of X 
meshes. 



1/3 



correspond to X x X x X fc-point 



B. Insulators 

We begin our computational study with a sequence of 
semiconductors and insulators in the high-symmetry di- 
amond structure. FIG. [2] compares the deviation of the 
exchange energy of silicon, cubic silicon carbide and di- 
amond at various finite fc-point configurations from the 
infinite limit for different singularity regularization meth- 
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FIG. 2. Convergence of exact and screened exchange energies 
for three semiconductors and insulators in the diamond (zinc- 
blende) structure: (a) silicon, (b) cubic silicon carbide (phase 
3C), and (c) diamond. See last paragraph of Section III A 



for details. The non-monotonicity in the absolute asymptotic 
error of the auxiliary-function results is due to a change in 
sign of that error. 



ods. 

While the order of magnitude of error in the exact ex- 
change energy with coarse fc-point meshes is comparable 
for all methods, Wigner-Seitz truncation typically yields 
significantly lower errors than do the other methods. 
The Wigner-Seitz truncated and the probe-charge Ewald 
methods (red +'s and blue *'s respectively in FIG. 2]) ex- 
hibit smooth convergence for all fc-point meshes including 
anisotropic ones, whereas the remaining methods incur 
higher errors for anisotropic fc-point meshes. The pat- 
tern of errors with fc-points for each method is similar 
for the three materials with the same underlying Bravais 
lattice and point-group symmetries. 

In contrast, the asymptotic exponential convergence 
of the truncated methods leads to orders of magnitude 
reduction in error for fine fc-point meshes, compared to 
the probe-charge Ewald and auxiliary-function methods, 
which exhibit L~ 3 convergence. The exponential decay 
length of the error in the exchange energy with respect to 
L, taken here to be the nearest neighbor distance in the 
effective fc-point sampled super-lattice, decreases from 
5.5 A in silicon through 3.5 A in cubic silicon carbide to 
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FIG. 3. Convergence of exact and screened exchange energies 
for hexagonal silicon carbide (phase 2H). See last paragraph 
of Section [HI Al for details. 
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FIG. 4. Convergence of exact and screened exchange ener- 
gies for proton-ordered cubic ice. See last paragraph of Sec- 
tion [TTT^] for details. 



2.5 A in diamond. The corresponding band gaps A are 
1.1 eV, 2.3 eV and 5.5 eV respectively. The decay length 
varies roughly as A -1 / 2 , similar to the density-matrix lo- 
calization length scale of tight-binding insulators.^ Con- 
sequently, the relative accuracy of the truncated meth- 
ods for insulators increases dramatically with increasing 
band gap as seen in FIG. [2j Note that the accuracy of 
the truncated potential methods for the exact exchange 
energy matches that of the screened exchange energy 
(red +'s and green x's versus pink A's), indicating that 
these methods are truly limited only by Brillouin-zone 
discretization errors in the underlying Kohn-Sham or- 
bit als. 

In high-symmetry materials, the accuracy of spheri- 
cal truncation is similar to Wigner-Seitz truncation for 
isotropic fc-point meshes since the Wigner-Seitz cell of the 
effective super-lattice is approximately spherical (rhom- 
bic dodecahedron for the FCC unit cell in the zinc-blende 
structure). This is no longer the case for lower sym- 
metry crystals such as hexagonal silicon carbide (phase 
2H with the wurtzite structure) shown in FIG. [3] In 
this anisotropic case, the accuracy of the Wigner-Seitz 
truncated method (red +'s) continues to match that 
of screened exchange (pink A's). On the other hand, 
spherical truncation (green x's), although still exponen- 
tially convergent, is an order of magnitude less accu- 
rate. Similarly, amongst the asymptotically L~ 3 conver- 
gent methods, the superior accuracy of the probe-charge 
Ewald method for anisotropic /c-point meshes in the high- 
symmetry crystal carries forward to superior accuracy 
overall for lower-symmetry crystals, in comparison to the 
auxiliary-function method. 

The differences between the methods are most dra- 
matic for proton-ordered cubic ice XIc (the proposed 
ground state structure^) shown in FIG. ffl The highly- 
localized states in this material cause dramatic improve- 
ments in accuracy for the truncated methods even for 
coarse fc-point meshes. Once again, lowered symmetry 
significantly favors the probe-charge Ewald method in 
comparison to the auxiliary- function method (blue ^-'s 
versus cyan CPs). 
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FIG. 5. Convergence of exact and screened exchange energies 
for graphite. See last paragraph of Section [HI A| for details. 



C. Metals 

As demonstrated above, the exponential localization 
of Wannier functions leads directly to exponential con- 
vergence in the case of truncated Coulomb interaction 
methods. In contrast, we expect the discontinuity at the 
Fermi surface at zero temperature to lead to algebraic 
convergence in metallic systems, which we explore now. 

FIG. [5] shows the convergence behavior of the various 
methods for the case of graphite, which is semi-metallic. 
Once again, the probe-charge Ewald and auxiliary- 
function methods exhibit L~ 3 convergence. Wigner-Seitz 
truncation no longer exhibits exponential convergence, 
but remains the most accurate method for computing 
the exchange energy, with accuracy comparable to that 
of screened exchange as before. Spherical truncation and 
the auxiliary-function method are less accurate due to 
the lower symmetry of the crystal structure in this case. 

The exponents governing the localization in graphite 
are complicated by the layered quasi-two-dimensional 
structure with weak inter-planar coupling. We ana- 
lyze those details in the related two-dimensional mate- 
rial graphene in Section |III D[ and here now focus on a 
simpler, three-dimensional metal, platinum. 

In simple metals, the Wannier-like functions wf {r) 
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FIG. 6. Convergence of exact and screened exchange energies 
for face-centered cubic metallic platinum at electron temper- 
atures, (a) T = 0.1 eV and (b) T = 1 eV. See last para- 
graph of Section |III A| for details. The non-monotonicity in 
the exponentially-convergent results in (b) for exact exchange 
with truncated-potentials and screened exchange is due to a 

1 /3 

change in sign of the error near N k ~ 12. 



FIG. 7. Convergence of exact and screened exchange energies 
for graphene at electron temperatures, (a) T = and (b) 
T = 1 eV. See last paragraph of Section [HI A| for details. The 
non-monotonicity in the exponentially-convergent results in 
(b) for Wigner-Seitz truncated exact exchange is due to a 



1/2 

change in sign of the error near NJ ~ 11 



given by ^ decay 



as discussed in Section II A 



This leads to r decay of the pair densities and con- 



sequently ~ 4nr 2 dr(l/r)r 



L~ 2 errors due to 
truncation in the Coulomb self energies in (10). The 
L~ 2 errors dominate the asymptotic convergence of all 
the methods for metals at low temperatures, as shown in 
FIG. pla) for platinum. However, Wigner-Seitz trunca- 
tion (red +'s) continues to yield the highest accuracy for 
exact exchange in practice, particularly for coarse fc-point 
meshes. 

At finite Fermi temperature T for the electrons, the 
exponential decay length scales as at/T, where a is the 
lattice constant and t is the typical band width. When 
the number of A-points along each dimension exceeds ap- 
proximately t/T, this decay length plays an analogous 
role to the Wannier-function length scale of insulators. 
FIG. [6^b) shows the restored exponential convergence of 
the truncated methods and screened exchange, and the 
L~ 3 convergence of the other methods, in this regime. 



D. Lower dimensional materials 

In lower dimensional semiconducting or insulating sys- 
tems, we should expect the localization of the underlying 
Wannier functions to allow for exponential fc-point con- 
vergence with an appropriately chosen method from Sec- 



ality can lead to different exponents for the polynomial 
convergence, which we also explore here. 

The semi-metallic behavior of the two-dimensional ma- 
terial graphene is particularly interesting. In this case, 
the localization properties are determined by the phase 
twist of Bloch functions in fc-space about the Dirac point. 
The two-dimensional Fourier transform of that phase 
twist yields an r~ 2 decay of the Wannier-like functions. 
The pair densities fall off as r~ A leading to truncation 



tion II B For the metallic cases, the reduced dimension- 



errors in the Coulomb self energies in ( 10 ) that scale 
as 2nrdr(l/r)r- 4 - L~ 3 . FIG. [^a) shows that all 
methods, therefore, exhibit L~ 3 asymptotic convergence 
for graphene at zero temperature. The errors oscillate 
with a period of 3 fc-points per dimension because the 
discrete fc-point mesh includes the special Dirac point 
when the sampling is a multiple of 3. 

At sufficiently high temperatures, the exponential 
length scale at/T of the Wannier-like functions becomes 
relevant at practical fc-point meshes, just as in three- 
dimensional metals. FIG. shows that this length 
scale restores the exponential convergence of the Wigner- 
Seitz truncated method and screened exchange (red +'s 
and pink A's respectively). 

Truncation on a three-dimensional sphere is no longer 
meaningful in these lower dimensional materials, and 
it gives reasonable results only for intermediate fc-point 
meshes which minimize the aspect ratio of the super- 
cell. Analytically Fourier transforming the Coulomb po- 
tential truncated on the appropriate 'lower-dimensional 
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FIG. 8. Convergence of exact and screened exchange energies 
for the semiconducting (8,0) single-walled carbon nanotube. 
See last paragraph of Section [HI A| for details. 



spheres', finite cylinders in two-dimensional materials 
and finite right prisms in one-dimensional materials, is 
no longer possible. Wigner-Seitz truncation using the 



MIC algorithm (A4) is clearly the method of choice in 
these geometries. 

Finally, in one dimensional systems, the probe-charge 
Ewald and auxiliary-function methods are identical as 
proved in Section [II B| As FIG. [8] shows, for a semicon- 
ducting (8,0) single- walled carbon nanotube (SWCNT), 
both of these methods yield rather poor L~ 3 convergence, 
in contrast to the exponential convergence of the Wigner- 
Seitz truncation and the screened exchange interaction 
for this system. 



E. Total Energy Convergence 

The results in the preceding sections demonstrate that 
calculation of exact exchange with Wigner-Seitz trunca- 
tion, with relatively few exceptions, generally requires 
fewer fc-points to reach a given level of convergence than 
all other methods. Moreover, we have seen that cal- 
culation of the long-ranged exact exchange, when per- 
formed with the Wigner-Seitz truncated method, com- 
petes with the short-ranged screened exchange of the 
HSE06 functional, which models the long-ranged compo- 
nents of the exchange energy within a semi-local approx- 
imation. We now ask whether the Wigner-Seitz method 
makes it possible to evaluate exact-exchange functionals 
on the same, relatively modest fc-point meshes needed for 
simple density-functional theory calculations. 

To address the above issues, FIG. [9] compares the total 
energy convergence of a purely semi-local density func- 
tional (PBFj23), a standard hybrid functional (PBEC 7 ) 
employing exact exchange computed using various stan- 
dard approaches as well as our Wigner-Seitz approach, 
and a hybrid functional (HSEOfP^) employing short- 
ranged screened exchange. The results in FIG. [9] show 
that the total energy convergence of the exact-exchange 
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FIG. 9. Total energy convergence of the semi-local functional 
PBE and the screened-exchange hybrid functional HSE06 
compared to the hybrid functional PBE0 with exact exchange 
computed using various methods for (a) silicon (b) hexago- 
nal silicon carbide (phase 2H) and metallic platinum at (c) 
T = 0.1 eV and (d) T = 1 eV. The logarithmic energy scale 
shows the deviation of the energy for each functional at finite 
fc-point meshes from the converged energy for the same func- 
tional. The non-monotonicity in the exponentially-convergent 
results in (d) for hybrid functionals with truncated-potential 
exact exchange and screened exchange is due to a change in 

1/3 

sign of the error near N k ~ 12. 



hybrid functional PBE0, when computed using truncated 
potentials, is indeed comparable to that of the screened- 
exchange hybrid functional HSE06. Moreover, comput- 
ing exact exchange with truncated methods not only 
matches the convergence of traditional semi-local den- 
sity functionals, but sometimes even outperforms their 
convergence for insulators (as in FIG.[9£i,b). 

Spherical truncation yields similar convergence to 
Wigner-Seitz truncation for high-symmetry insulators 
(green x's versus red +'s in FIG(9k), and is less accu- 
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rate for lower symmetry ones (FIG. |9p), as expected. 
The auxiliary function and probe-charge Ewald methods 
limit the total energy convergence of PBEO to L~ 3 (cyan 
□ 's and blue sfc's in FIG. [9|a,b,d)), in contrast to the 
exponential convergence of PBE and HSE06. 

The r~ 2 localization of the Wannier-like functions in 
metals at low temperature limit the convergence of all 
methods, including the semi-local functionals (FIG. [9j:). 
The L~ 2 errors in the tails of the Kohn-Sham orbitals 
lead to L~ 2 errors in the Kohn-Sham kinetic energy, 
which dominate in this situation. Consequently, all meth- 
ods for treating exchange yield total-energy convergence 
of PBEO similar to that of PBE and HSE06 for low- 
temperature metals, with Wigner-Seitz truncation only 
marginally better than the others. 

At sufficiently high temperatures (FIG. |9ji), the ex- 
ponential length scale at/T becomes relevant at practi- 
cal fc-point configurations as before, leading to exponen- 
tial total-energy convergence of the semi-local function- 
als. As in the case of insulators, the auxiliary-function 
methods limit the accuracy of the hybrid functional while 
the truncated-potential methods yield total-energy con- 
vergence for the exact-exchange hybrid functional on par 
with the screened-exchange and traditional semi-local 
functionals. 



Coulomb potentials truncated on a subset of lattice di- 
rections. 

Finally, we explore the accuracy of several meth- 
ods for computing exchange energies applied to mate- 
rials with varying electronic structure and dimension- 
ality. The Wigner-Seitz truncation motivated in this 
work systematically yields the most accurate results for 
all of these systems, demonstrates exponential conver- 
gence with Brillouin-zone sampling for all but metals 
at zero temperature, and, most importantly, delivers k- 
point convergence for hybrid functionals on par with 
that of screened-exchange and even traditional semi-local 
functionals. The findings of this work will enable accu- 
rate calculations of periodic systems employing exact- 
exchange hybrid functionals at the same computational 
effort as those employing screened-exchange functionals, 
and bring us one step closer to widespread ab initio stud- 
ies of processes such as catalysis at solid surfaces that 
require chemical accuracy. 

We thank K. A. Schwarz and D. Gunceler for their 
helpful suggestions in improving this manuscript. This 
work was supported as a part of the Energy Materials 
Center at Cornell (EMC 2 ), an Energy Frontier Research 
Center funded by the U.S. Department of Energy, Office 
of Science, Office of Basic Energy Sciences under Award 
Number Ide-scOOO 10861 



IV. CONCLUSIONS 

Hybrid density-functionals enable high-accuracy pre- 
dictions within Kohn-Sham theory, but state-of-the-art 
methods for computing such functionals necessitate dense 
fc-point meshes in calculations of periodic systems due 
to the regularization methods commonly employed for 
the singular Coulomb integrals in the exact-exchange en- 
ergy. In this work, we analyze the exchange energy in 
a real-space formalism based on Wannier functions and 
show that the dominant errors at finite fc-point meshes 
arise from the resultant artificial periodicity of the Wan- 
nier functions. Truncating the Coulomb potential on the 
Wigner-Seitz cell of the effectively-sampled super-lattice 
is therefore the ideal method to minimize these errors. 

We then prove analytically that the exchange energy 
computed using truncated potentials converges exponen- 
tially with the number of fc-points, N^, whenever the one- 
particle density matrix is exponentially localized. This 
is the case for all but metallic systems at T = 0. In 
contrast, the frequently-employed auxiliary-function reg- 
ularization methods for the exchange kernel suffer from 
a dipole-dipole interaction with the artificial images that 
limits the asymptotic convergence to NZ . 

To deliver this exponential asymptotic convergence to 
practical calculations, we then develop in the appendix, 
general and efficient constructions for the Wigner-Seitz 
truncated Coulomb potentials in the plane-wave ba- 
sis. Additionally, we generalize the truncated-potential 
and auxiliary-function methods to slab-like or wire- 
like systems with lower-dimensional periodicity by using 



Appendix A: Truncated Coulomb potentials 

Computations for periodic systems frequently employ 
the plane-wave basis™ which presents the advantage of 
systematic and exponential basis-set convergence con- 
trolled by a single parameter, namely the kinetic energy 
cutoff, or equivalently, the Nyquist frequency. This ad- 
vantage can be extended to non-periodic systems and sys- 
tems with lower-dimensional periodicity such as slabs and 
wires by using truncated Coulomb potentials! 35 * 36 ! Ad- 
ditionally, above, truncated Coulomb potentials proved 
particularly useful in the computation of the exchange 
energy even for periodic systems. 

As shown above, truncation of the potential on the 
Wigner-Seitz cell leads to the most accurate method for 
the exchange energy, and it is also the most efficient 
method for lower-dimensional geometries since, as dis- 
cussed by Ismail-BeigiJ^ it localizes the G = singu- 
larity of the Coulomb kernel to a single point. However, 
the singular Fourier integrals required to construct plane- 
wave Wigner-Seitz truncated Coulomb kernels cannot be 
solved analytically in general, and are prohibitively ex- 
pensive to compute numerically. Here, we develop a 
general and efficient O(NlnN) construction for these 
kernels based on the minimum-image convention (MIC) 
methodP^l Finally, as a practical matter, we note that 
determining the Coulomb kernel is a one-time computa- 
tion because the corresponding storage requirements are 
modest. 

To establish our normalization conventions for the 
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plane-wave basis, we expand periodic charge densities as 

pif) = J2g el< ~ ' r P<3 wnere G are reciprocal lattice vectors 
for the unit cell of volume Q. The interaction energy un- 
der a translationally invariant potential 



/':_» := J dn I df 2 p* 1 (r 1 )K(T f 1 ~T f 2)p2(r2) 



^HKg^g^G 



(Al) 



Employing a range-separation parameter a, we ap- 
proximate this kernel by 



kT = I df c -^( edcar i crfQr 

\ r r 

-g 2 \ n 



ws 

Si 1 exp 1 V, • 



rEWS 



=s j?erfar 



-iG-r 



(A4) 



is then diagonal in reciprocal space, with Kg = 

J dfe~ lGr K(r). Here, we denote integrals over unit cells 
by df and integrals over all space by J dr. For the 
long-ranged Coulomb interaction K(r) — 1/r, the plane- 
wave kernel Kg = Att/G 2 is singular at G = 0, and the 
interaction energy is finite only for neutral unit cells. 
In practice, the interaction energies of charged subsys- 
tems are computed by excluding the G = term, which 
amounts to adding a uniform neutralizing background 
charge to each subsystem. 

Treating non-periodic systems in the plane-wave basis 
requires the elimination of interactions between different 
unit cells in a translationally-invariant manner in order to 
preserve the efficiency and accuracy of the Fourier spec- 
tral method. If the Coulomb potential is truncated so 
that it is zero outside the first Wigner-Seitz cell and the 
charge densities are confined to half the domain of trun- 
cation, then the interaction within each unit cell remains 
unmodified and the interaction between unit cells is ex- 
actly zero. The simplest choice for accomplishing this 
truncates the Coulomb potential outside a sphere of ra- 
dius R smaller than the in-radius of the Wigner-Seitz cell, 
leading to the analytical plane-wave kernel, 



K s ~ h = 

G 



dre 



-iG-r^ 



(2ttR 2 



\r\<R. 



^(1 - cos GR), 



G = 
G> 0. 
(A2) 



where the discrete sum over f is a quadrature on the 
Wigner-Seitz cell with Np nodes as described below. The 
short-ranged first term is localized to the Wigner-Seitz 
cell by choice of a, so that it is unaffected by the trunca- 
tion and can be evaluated analytically. The error in this 
term can be reduced to machine precision e by choos- 
ing a = y — In e/i?i n , where R[ n is the in-radius of the 
Wigner-Seitz cell. The G = component is well-defined 
due to the finite real-space range, and is understood to 
be equal to its G — > limit given by ir/a 2 . 



The second term of ( A4 ) has a long-ranged smooth in- 



tegrand and is approximated by a Gauss-Fourier quadra- 
ture, evaluated as a fast Fourier transform (FFT) in prac- 
tice. This quadrature consists of nodes on a uniform par- 
allelepiped mesh with uniform weights, which we remap 
using the periodicity of the lattice to the first Wigner- 
Seitz cell. In the interior of the integration domain, the 
integrand is bandwidth-limited as exp(— G 2 /4a 2 ), and 
the error in the Fourier quadrature can be reduced to 
e by choosing an FFT resolution such that the Nyquist 
frequency exceeds 2ay — lne = — 21ne/i?j n . The cusps 
in the periodic repetition of the integrand at the bound- 
aries of the Wigner-Seitz cell cause an additional error 
in the kernel, but this error does not contribute in the 
Coulomb energy of charge distributions that are confined 
to the half-sized Wigner-Seitz cell and are resolvable on 
the chosen Fourier gridP2l 



1. Minimum-image convention (MIC) method 

From the above arguments, the natural choice for the 
truncated potential is clearly 

Kj s = [ dfe- iSf - } (A3) 
G Jws r 

where J wg represents integration over the first Wigner- 
Seitz cell. However, the singularity at the origin and the 
polyhedral domain of integration preclude general ana- 
lytical solutions and straightforward numerical quadra- 
tures. Martyna and Tuckerman introduced an approxi- 
mate construction for Coulomb potentials truncated on 
parallelepiped domains based on range-separation tech- 
niques. Here, we generalize this so-called minimum- 
image convention (MIC) methocP^ to Wigner-Seitz cells 
of arbitrary lattice systems. 



2. Partially-truncated Coulomb kernels 

Next, we generalize the above construction to systems 
with lower-dimensional periodicity, where the Coulomb 
kernel is truncated along some lattice directions and re- 
mains long-ranged along the others. In these geometries, 
the kernel is still singular around G = 0, albeit with a 
slower divergence: In G for one periodic direction or 1 jG 
for two periodic directions in contrast to 1/G 2 for the 
fully periodic case. A general shape for the truncation 
domain in the non-periodic directions leads to a kernel 
which is singular for an entire line or plane of recipro- 
cal lattice vectors passing through G = 0. Ismail-BeigP^ 
pointed out that Wigner-Seitz truncation localizes the 
singularity to the single point G = in all these cases. 

We can understand the special property of the Wigner- 
Seitz truncation by writing the Coulomb kernel truncated 
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on an arbitrary domain D as 



/V;' = / dre- l6 ?l -B D (r) = f dk^6 D (G-k), (A5) 



G 



for any two-dimensional lattice systemP^ Accordingly, 
we generalize the MIC approach^ anc j approximate the 
partially-truncated wire-geometry Coulomb kernel by 



T>-wirc 
G 



47T 

G 2 



1 — exp 



-G 2 
Aa 2 

r± ?i6WSj 



-iG-r± (~iol 



Cfc m{ (r ± ) (A8) 



where r± are nodes for the two dimensional Gauss- 
Fourier quadrature mapped down to the Wigner Seitz 
cell, WSj_ , of the truncated directions with area £l± . 

Here, we introduce the smooth, long-ranged special 
function C%(p), which plays the same role for Ck(p) 
that eii(ar)/r plays for 1/r in the fully-truncated case 
of (A4). Operationally, this function is defined by the 



two-dimensional convolution 



where Or)(f) is a function that is 1 for r € D and 
otherwise, and Od is its Fourier transform. Since Ooif) 
is constant along periodic directions, 9 is zero for wave- 
vectors with any component along those directions. In 
general, the singularity from 4n/k 2 'infects' all wave- 
vectors with no component along the periodic directions. 
Hence, the singularity is spread to a plane of points for 
one-dimensional or wire-like systems, and a line of points 
for two-dimensional or slab-like systems. However, the 
Fourier transform of a ^-function with a shape that tiles 
with the periodicity of the lattice, such as the Wigner- 
Seitz cell, is zero at all non-zero reciprocal lattice vectors. 
This confines the singularity to G — for these special 
cases. 

Now consider, without loss of generality, the slab ge- 
ometry with its one truncated direction along z, and let L 
be the unit cell length along that direction. The Fourier 
transform of 1/r over the two periodic directions eval- 
uates to 2ne~ Gp \ z \ /G p where G p is the component of 
the wave-vector along the untruncated directions. At 

G p = 0, removing the singular part 2n/G p leaves be- For k = 0, (A9| reduces to the analytical^ expression 
hind — 27r|z|, the potential due to an infinite plane of 
charge with the arbitrary offset in potential fixed to be 
zero at the plane. Note that, although this choice of 
zero of potential does not change the total energy for 
a neutral charge distribution, care must be exercised to 
ensure that it be consistent for all interactions between 
charged subsystems of an overall neutral system. Finally, 
the remaining integral over the Wigner-Seitz cell in the 
truncated direction, i.e. zG [-L/2, L/2), can also be 
performed analytically) 35 * 36 ! so that the Coulomb kernel 
for truncation in slab geometries becomes 



2a 2 p , dp'e- a2 (p 2 +^I (2a 2 pp>)C k (p'). 

(A9) 



C$(p) = -21np - T (a 2 p 2 ), but for k ^ 0, C£(p) needs 
to be parametrized numerically.^ The choice of a and 
FFT resolution in ( A9 ) follow the discussion for the fully- 
truncated MIC construction, except that R in is the radius 
of the two dimensional Wigner-Seitz cell WS_l , and inde- 
pendent two-dimensional fast Fourier transforms produce 
the results for each plane of constant G z . 



G 



lab 




G Z L 



exp ■ 



-G„L 



G^Q 
G = 0. 



(A6) 



Similarly, for the wire geometry with the single pe- 
riodic direction along z, the partial Fourier transform 
of 1/r over the periodic direction is 2Ko(G z p), where 
p = v 'r 2 — z 2 is the usual cylindrical coordinate, and 
Kq is the modified Bcsscl function of the second kind. 
At G z = 0, removing the logarithmically divergent part 
leaves behind —2 In p so that the regularized partial 
Fourier transform of the Coulomb potential at G z = k is 



3. Ewald sums for reduced-dimensional systems 

The plane-wave Coulomb kernels above, truncated over 
the Wigner-Seitz cell in one, two or three lattice direc- 
tions, enable the calculation of Coulomb interaction en- 
ergies in slab, wire and isolated geometries respectively. 
However, a purely reciprocal-space method is only prac- 
tical if at least one of the two charge densities, pi(r) 
or p2(r) in (Al), is bandwidth limited. The interaction 
energy of point nuclei with each other does not satisfy 
this criterion and requires the use of an Ewald sumPSl 
Generalizing the standard Ewald method to an arbitrary 
combination of truncated and periodic lattice directions, 
gives the interaction energy for a set of point charges Zi 
at locations r,- in the first unit cell, as 



2A 0(M , k*0 
yH> l-21np, fc = 0. 



(A7) 



E, 



ewald 



The remaining Fourier transform over the two truncated 
directions is analytically computable for a cylindrical 
truncation domain,^ but that choice spreads the diver- 
gence beyond G = as mentioned previously. On the 
other hand, the Fourier transform of ( |A7[ ) with a Wigner- 
Seitz truncation domain is not known in closed form 



E 



R,i,j 
i^Lj if R=0 

z,z 



ZiZj erfc r)\fi+R- 
2 \n + R-fA 



y ^ iG . ( ._ f , ) _ _ny Z 2 



G,i,j 



per 



(A10) 
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Here, the first term evaluates the contribution due to the 
short-ranged part erfc(?yr)/r of the Coulomb potential, 
the second term captures the contribution due to the re- 
maining long-ranged part erf(r}r)/r, and the third term 
exactly cancels the self interactions introduced by the sec- 
ond term. The standard range-separation parameter r\ is 
adjusted to simultaneously optimize the convergence of 
the sum over lattice vectors R as well as that over recip- 
rocal lattice vectors G. (When some lattice directions are 
truncated, R and G correspond to the lattice vectors and 
reciprocal lattice vectors of the lower dimensional Bravais 
lattice of periodic directions alone.) Finally, in the second 
term, f2 pcr is the volume, area or length of the unit cell 
along the periodic directions alone and g%{r) is related 
to the (partial) Fourier transform over those directions 
of the long-ranged part of the Coulomb potential. 



When all three directions are periodic, g v -(r) 



-G 
4r) 2 



the double sum over point charges factorizes 



exp 

to the square of the structure factor, and (|A10|) reduces 



to the standard Ewald sumPSl Next, for the slab geom- 
etry truncated, without loss of generality, along the z 
direction, we find 



4(r) 



— 2tt (^z crf(?7z) + 



G^O 
G = 0, 



(All) 



where fcj{z) = e Gz erfc (G/2r] + rjz), and this reduces 



(A10I to the 'Ewald 2D' formula 



38 39 



The Ewald sum for the wire geometry with one peri- 
odic and two truncated directions does not seem to have 
been addressed previously, perhaps because g^ is not an- 
alytically expressible in that case. In fact, we can show 
that gl(r) = C G (Vr> 



z 2 ), precisely the function de- 
fined in (A9|, which was introduced for our generaliza- 



tion of the MIC method to this geometry. Finally, when 
all three lattice directions are truncated, the Coulomb 
kernel has no G = singularity, and an Ewald sum is 
not required. In this case, the Coulomb energy of a set 
of point charges is computed directly in real space as a 
sum over all pairs in one unit cell. 
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